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Advanced numerical methods in astrophysical fluid 

dynamics 



Computational gas dynamics has become a prominent research field both in as- 
trophysics and cosmology. In the first part of this review we intend to briefly 
describe several of the numerical methods used in this field, discuss their range 
of application and present strategies for converting conditionally-stable numerical 
methods into unconditionally-stable solution procedures. The underlying aim of 
the conversion is to enhance the robustness and unification of numerical methods 
and subsequently enlarge their range of applications considerably. In the second 
part Fabian Heitsch presents and discusses the implementation of a time-explicit 
MHD Boltzmann solver. 



1.1 Numerical methods in AFD 

Astrophysical fluid dynamics (AFD) deals with the properties of gaseous-matter 
under a wide variety of circumstances. Most astrophysical fluid flows evolve over a 
large variety of different time and length scales, henceforth making their analytical 
treatment unfeasible. 

On the other hand, numerical treatments by means of computer codes has wit- 
nessed an exponential growth during the last two decades due to the rapid devel- 
opment of hardware technology. Nowadays, the vast majority of numerical codes 
are capable of treating large and sophisticated multi-scale fluid problems with high 
resolutions and even in three-dimensions. 

The numerical methods employed in AFD can be classified into two categories: 

(i) Microscopic oriented methods mostly based on N-body (NB), Monte-Carlo 
(MC) and on the Smoothed Particle Hydrodynamics (SPH). 

(ii) Grid oriented methods. To this category belong the finite difference (FDM), 
finite volume (FVM) and finite-element methods (FEM). 

Most numerical methods used in AFD are conditionally-stable. Hence, they may 
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converge if the Courant-Friedlichs-Levy condition for stability is fulfilled. As long 
as efficiency is concerned, these methods are unrivalled candidates for flows that 
are strongly time-dependent and compressible. They may stagnate however, if 
important physical effects are to be considered or even if the flow is weakly in- 
compressible. On the other hand, only a small number of the numerical methods 
employed in AFD are unconditionally stable. These are implicit methods, but they 
are effort-demanding from the programming point of view. 

It has been shown that strongly implicit (henceforth IM) and exp Ucit (hencefort h 
EM) methods aie different variants of the same algebraic problem (lHujeirati . l2005l) . 
Hence both methods can be unified within the context of the hierarchical solution 
scenario (henceforth HSS, see Fig. 11.31 ). 

In Table 11.11 we have summarized the relevant properties of several numerical 
methods available. 




Microscopic oriented 
methods 



NB, MCM, SPH 



Conditionally stable] 
Methods 



(Unconditionally stable 
V Methods 




%. R< Jch<_< Vf<_< 5hp<< _"TtS<<_" 5yii<< 



'Acq 



Fig. 1.1. Numerical methods (finite difference, finite volume, finite element, N-Body, 
Monte Carlo and the smoothed particle hydrodynamics employed in AFD and their possi- 
ble regime of application from the time scale point of view. The time scales read as follows: 
the radiative-TR, gravitative-Tc, chemical-Tch, magnetic-TMp, hydrodynamic-THD, 
thermal-Txh, viscous-ryis, and the accretion time scale-TAcc- 



1.2 Time scales in AFD 

Assume we are given a box of Lx Lx L dimensions filled with a rotating multi- 
component gaseous-matter. The gas is said to be radiating, magnetized, chemical- 
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Explicit Implicit HSS 



solution 
method 



^"+1 = q" +6td" 



q"+^ ^q" +5tA-^d* 



: m/' + (1 -a)6tA-/d' 



Type of flows 



Strongly time- 
dependent, 
compressible, 
weakly dissipative 
HD and MHD 
in 1, 2 and 3 dimen- 
sions 



Stationary, 
quasi-stationary, 
highly dissipative, 
radiative and 
axi-symmetric 
MHD-flows in 1, 
and 3 dimensions 



Stationary, 
quasi-stationary, 
weakly compressible, 
highly dissipative, 
radiative and 
axi-symmetric MHD- 
flows in 1, 2 and 3 di- 
mensions 



Stability 


conditioned 


unconditioned 


unconditioned 


Efficiency 


1 (normalized/2D) 


2 

~ m 


2 

~ m, 

d 


Efficiency: 

Enhancement 

strategies 


Parallelization 


Parallelization, 

preconditioning, 

multigrid 


HSS, parallelization, 

preconditioning, 

prolongation 


Robustness: 

Enhancement 

strategies 


i. subtime-stepping 

ii. stiff terms 
are solved 
semi-implicitly 


i. multiple iteration 

ii. reducing the time 
step size 


i. multiple iteration 

ii. reducing the time 
step size, HSS 


Numerical 

Codes 

Newtonian 


Solvers 1" 

ZEUS&ATHENA*, 
FLASHY NIRVANA'', 
PLUTO^ VAC ^ 


Solver2« 


IRMHD'' 


Numerical 

Codes 

Relativistic 


Solvers3' 

GRMHD^, ENZO*, 
PLUTO', HARM'", 
RAISHIN",RAM°, 
GENESIS'', WHISKY 


Solver4'' 


GR-I-RMHD' 



Table 1.1. A list of only a part of the grid-oriented codes in AFD and their 
algorithmic properties. In these equations, q"'"'^^, 6t, A, a and d* denote the 

vector of variables from the old and new time levels, time step size, a 
preconditioning matrix, a switch on/off parameter and a time-modified defect 
vector, respectively, "m " in row 4 denotes the bandwidth of the corresponding 

matrix. 

~ Clarke (1996), ' Stone, Norman ('1992'); 



'^Bodenheimer et al. (1978); Clarke (1996), ' Stone, Norman (1992); Gardiner, Stone (2006). 'Frvxel l et alJ 
(2000), ' Ziegler cl998.) , 'Mignone Bodo (2003 ); Mignone et al. (2007), 'Tothetal. (1998), 'WuchteJ 

t]); Swestv <1995t). 'iHuieii-all {\993l |2005|); i Falle (2003), Koide et al. (1999); Komissarov (20MX 
Villiers, HawlevI 120031), 10'Sheaet alJ )2"o"04.), ■.Mignone et al., (.2007.), ".Gammie et aU t2003,) 
' Mizuno et al. ( 2006), ' Zha ng. MacFadveij <2006l) . lAlav et alj<1999l) . iBaiotti et alji2003ft , iLiebendorfer et alj 
gOOZ) . ' Huieiratet al. . ( ,20071) ! 
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reacting, partially ionized and under the influence of its own/external gravitational 
tield. Let the initial state of the gas be characterized by a constant velocity, density, 
temperature and a constant magnetic field. The time-scales associated with the 
flow can be ob tained directly from the radiative MHD-equations as follows (see 



Hujeiratll2005L for detailed description of the set of equations). 




Fig. 1.2. The regime of application of explicit method is severely limited to Euler-type 
flows, whereas sophisticated treatment of most flow-problems in AFD require the employ- 
ment of much more robust methods. 



• Continuity equation: 



^ + V-pV = 0, (1.1) 
ot 



where p, V stand for the density and the velocity field. Using scaling variables 
(e.g. Table ll.2l ). we may approximate the terms of this equation as follows: 

— ~ — and V-pV ~ This yields the hydrodynamical time scale thd - —■ 
dt T L V 

The so-called accretion time scale can be obtained by integrating the continuity 
equation over the whole fluid volume. Specifically, 

r —dVol = — ~—, f {V-pV)dVol= f pV-n-dS ^ AM ~ M, 

Jvol Ot dt T Jvol Js 
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Scaling variables Molecular cloud Accretion(onto SMBH) Accretion (onto UCO) 



I 


Length 


0(pc) 


<9(AU) 


(9(10^- cm) 


p 


Density 


10^22 g cm--' 


lO"*" g cm-^ 


10"^ gcm"^ 


T 


Temperature 


10 K 


10^ K 


10^ K 


V 


Velocity 


0.3 kms"' 


10^ km s"' 


10^-^ km s"' 


B 


Magnetic Fields 


30 pG 


lO^G 


lO^G 


M 


Mass 


10^ Mo 


lO*- Mo 


Mo 


M 


Accretion rate 




10-2 Mo 


10-'° Mo 



Table 1.2. A list of possible scaling variables typical for three different 
astrophysical phenomena: giant molecular clouds, accretion onto supermassive 
black holes (SMBHs) and accretion onto ultra-compact objects (UCO). These 
variables may be used for reformulating the radiative MHD equations in 

non-dimensional form. 



where "Vol" denotes the total volume of the gas and "S" coiTesponds to its sur- 
face. Equating the latter two terms, we obtain: 

M . M 

M ^ Tacc ~ -r- 

T M 

In general Tacc is one of the longest time scales characterizing astrophysical 
flows connected to the accretion phenomena. 
The momentum equations: 



+ VV ® V - —VP + fccn, + ^ + V^ + — + Ql,^, (1.2) 

ot p p 4np 

where P, fcent, fad, ip, B, denote gas pressure, centrifugal force, radiative 
force, gravitational potential, magnetic field and viscous operators, respectively. 
From this equation, we may obtain the following time scales: 

(i) The sound speed crossing time can be obtained by comparing the follow- 
ing two terms: 

V ■ 1H 

— ~ — , which yields: r, ~ thd — 

dt p \Vs 

where Kis the sound speed. 

(ii) The gravitational time scale: 

where = GM/L and G is the gravitational constant. 



1.2 Time scales in AFD 
(iii) Similarly, the Alfven-wave crossing-time: 

dV VxBxB (V^ 



dt ' Anp ■ """^ \Va, 

where V^(= B^/Anp) denotes the Alfven speed squared. 

(iv) Radiative effects in moving flows propagate on the radiative scale, which 
is obtained from: 

dv u _ (yf 

-r- « ^ Trad - THD — , 

Ot p \c I 

where c is the speed of light. 

(v) The viscous time scale: 

where v is a viscosity coefficient. 

• The induction equation, taking into account the effects of o-dyn-dynamo, mag- 
netic diffusivity v^iff and of ambipolar diffusion reads: 

dB B 

— ^ V X {V X B + adynB - v^^gV X B) + V X I- x[Bx{VxB)]}, (1.3) 

dt 4nypiPn 

where pi^n denote the ion and neutral densities. 

Thus, the induction equation contains several important time scales: 
(i) The dynamo amplification time scale, which results from the equality: 
dB ^ ^ L 

— = V X adynS ^ Tdyn 



dt ttdyn 

(ii) The magnetic-diffusion time scale: 

dB L?' 

— ^ V X (VmagV X S) ^ Tdiff = 

Ot ^mag 

(iii) The ambipolar diffusion time scale: 

^ ^ V X { ^ X [B X (V X B)]] 

B 1 / BM/ 1 Vl B _ B 

- ~ T 1 T T^= ^amb-TT =^ T^mb = 



T L\Anpnj\ypij\LI ypi L'^ " Damb 
where Dambi- V'^/(yPi)) is the ambipolar- diffusion coefficient. 



10 



Advanced numerical methods in astrophysical fluid dynamics 



The chemical reaction equations. 

The equation describing the chemical-evolution of species "/" is : 



dt 



^ ^ hnnPrnPn + ^ ImPm , (1-4) 



where ^„,„ denotes the reaction rate between the species m and n. I,„ stands for 
other external sources. For example, the reaction equation of atomic hydrogen 
in a primordial gas reads: 

dpH k2 h ph k2 rriH 

-PH+Pe PHPe <^ PH Pe ^ Tch 



ot niH niH T mn K2 Pe 

where pe, ^2(10^° cm^ s~^) correspond to the electron density and to the gener- 
ation rate of atomic hydrogen through the capture of electrons by ionized atomic 
hydrogen, nin corresponds to the mass of atomic hydrogen. 
Equations of relativistic MHD 

The velocities in relativistic flows are comparable to the speed of light. This 
implies that the hydrodynamical thd and radiative T,-ad time scales are compa- 
rable and that both are much shorter than in Newtonian flows. 



Time scales Molecular cloud Accretion(onto SMBH) Accretion (onto UCO) 



Thd 


~ IQf' Yr 


~ months 


~ 1 s 


Trad/ Thd 


~ lO-*" 


~ 10-3 


~ 10-3 


Tgrav/THD 


~ 10-2 


~ 10-3 


~ 10-3 


Tch/thd 


~ 10-' 


~ 10-5 


~ 10-4 


T magi Thd 


~ 10-2 


~ 10" 


~ 10-' 


TvislTHD 


~ 10' 


~ 102 


~ 102 


TacclTHD 




~ 10^ 


~ 10'2 



Table 1.3. A list of the time scales relative to the hydrodynamical time scale for 
three different astrophysical phenomena. 



We note that although the dynamical time scale in relativistically moving flows 
is relatively short, there are still several reasons that justify the use of implicit 
numerical procedures. In particular: 

(i) The relativistic MHD equations are strongly non-linear, giving rise to fast 
growing non-linear perturbations, imposing thereby a further restriction on 
the size of the time step. 

(ii) The deformation of the geometry grows non-linearly when approaching 
the black hole. Thus, in order to capture flow-configurations in the vicinity 
of a black hole accurately, a non-linear distribution of the grid points is 
necessary, which, again, may destabilize explicit schemes. 



1.3 Numerical methods: a unification approach 
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(iii) Initially non-relativistic flows may become ultra-relativistic or vice versa. 
However, almost all non-relativistic astrophysical flows known to date are 
considered to be dissipative and diffusive. Therefore, in order to track their 
time-evolution reliably, the employed numerical solver should be capable 
of treating the corresponding second order viscous terms properly. 

(iv) The accumulated round off errors resulting from performing a large num- 
ber of time-extrapolations for time-advancing a numerical hydrodynamical 
solution may easily cause divergence. The constraining effects of boundary 
conditions may fail to configure the final numerical solution. 



Hierarchical solution scenario 




Fig. 1.3. A schematic description of the hierarchical solution scenario (HSS). The HSS 
is based on dynamical-varying the efficiency and robustness of the numerical method to 
leapfrog the transient phase. The method is most suitable for searching quasi-stationary 
flow-configurations that depend weakly on the initial conditions. Here the coupling be- 
tween the equations can be enhanced gradually, by starting solving them sequentially, then 
partial-coupling in combination with the operator splitting approach (OSA), full-coupling 
using the Krylov-subspace iterative method (KSIM) and finally extending the coupling to 
include the radiative transfer equation (RTE) and energy equation of multi-temperature 
plasmas. 



1.3 Numerical methods: a unification approach 

In this section we show that explicit and implicit methods are special cases of a 
more general solution method in higher dimensions. 
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Assume we are given the following evolution equation of a vector variable q: 

^+L{q) = f, (1.5) 
ot 

where L, / correspond to an advection operator and to external forces. 

Adopting a time-forward discretization procedure, the unknown vector q at the 
new time level can be extrapolated as follows: 



tf+' =tf' +6t-RHS", (1.6) 

where RHS = f-L{q). 

Depending on the time step size and on the number of grid points, the numerical 
procedure can be made sufficiently accurate in space and time. 
On the other hand Equation |l.6| can be viewed as an equality of two one-dimensional 
vectors: 

[vector of unknowns] = [vector of knowns] o q"^^ =b, (1.7) 
where ^ - q"+St-RHS". 

In higher dimensions, however. Equation 1 1.71 is a special case of the matrix equa- 
tion: 

Aq"^^ - b, (1.8) 

in which it is projected along the diagonal elements. It is obvious that the matrix 
I/6t is a further simplification of the matrix that contains just the diagonal elements 
of A. 

Therefore, we may adopt the higher dimension formulation to gain a better un- 
derstanding of the stability of the solution procedure. According to matrix algebra, 
a necessary condition for the matrix A to have a stable inversion procedure is that 
A must be strictly diagonally dominant. Equivalently, the entries in each row of 
the matrix A must fulfill the following condition: the module of the diagonal ele- 
ment dij is larger than the sum of all off-diagonal elements Xj^i where i and 
j denote the row and column numbers of the matrix . Applying a conservative 
and monotonicity preserving scheme, the latter inequality may be re-written in the 
following form: 

1 



— -I- positive contributions 



>XiM- (1-9) 

j(#i) 



We note that since 6t is a free parameter, it can be chosen sufficiently small, so 
that l/dt largely dominates all other off-diagonal elements, or so large that l/6t 
becomes negligibly small. 



1.3 Numerical methods: a unification approach 
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I ' ' ' ' ' ' ' ' ' 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

X 



Fig. 1.4. The the profile of the shock tube problem obtained with Courant-Friedrichs- 
Levy numbers CFL=0.4 and 0.9 using the PLUTO code. Although both CFL-numbers 
are smaller than unity the numerical solution procedure does not appear to be stable even 
with CFL=0.9. 

We may further simplify this inequality by choosing the time step size even 
smaller, such that 

^>J]kj|, forVj^/ (1.10) 

can be safely fulfilled. We may decompose the matrix A as follows: 

A = D + R = D{I + D'^R), 

where D is the matrix consisting of the diagonal entries of A and R {- A - D) 
consists of the off-diagonals. Thus, the elements of D are proportional to 1 16t, 
whereas those of R are proportional to 6t. This implies that A can be expanded 
around II6t in the form: 

A = A(°) + A^'^ + A^^) + • • •, (1.11) 

where the leading matrix A^^^ ~ I ^i^d ^^^^ ~ 5t I. In this case the inversion of 
the matrix A is not more necessary and the resulting numerical procedure would 
correspond to a classical time-explicit method. 
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1.3.1 Example 

The time-evolution of density in one-dimension is described by the continuity 
equation: 

The corresponding Jacobian matrix is: A = dLp/dp. The non-zero entries of A 
read: 



an- h and the off diagonal au - ~ — fori?!;] , (1-13) 

6t Ax Ax 



where Ax and /, j denote the grid spacing and grid point numbering. Applying 
a first order upwind discretization, then the condition of diagonal dominance de- 
mands: 



6t Ax 



(1.14) 



Ax 



This condition can be further simplified by choosing the time step size so small, 
such that 

1 2 max(|[/y|,|[/,+i|) 6tmax{\Uj\,\Uj^i\) 1 
— > o < -. (1.15) 

6t Ax Ax 2 

Thus, the condition of diagonal dominance is more restrictive than the normal 
CFL condition. This may explain, why most explicit methods fail to converge for 
Courant-Friedrichs-Levy number CFL = I - e (see Fig. 11.41 ). 




0.0 0.3 0.4 O.E O.B i.U 1.3 1<1[> lODO 1 DODO IDD IDOD 1 DOgO 

Radius Kum. af liwnB s^ops Mum. aF Mms s^Dps 



Fig. 1.5. Weakly incompressible flow between two concentric rotating spheres. Left panel: 
the 2D-distribution of Mach number (= V/V^) is displayed (25 isolines) for extreme 
weakly incompressible flows (Max (Mach) ~ 10"^). The maximum residual (middle) 
and the CFL-number (right) versus the number of time steps are shown. 



1.4 Converting time-explicit into implicit solution methods 
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1.4 Converting time-explicit into implicit solution methods 

In a series of publications, we have shown that the robustness of exphcit meth- 
ods can be enh anced gradually to recover full-implicit solution procedures (see 
Hujeirai 120051) . In the following we outline the main algorithmic steps towards 
extending classical explicit methods into implicit: 

(i) Use the same mathematical form of RHS" of Eg. |1.6| to compute RHS"'*'^ = 
RHS iq"+^) and subsequendy the mean RHS - a ■ RHS " + {l-a)- RHS 
where < a < 1 is a pai'ameter that may depend also on the time step size. 

(ii) Define the defect 



n+l _ n 

d = -{- —) + RHS. 

6t 



(1.16) 



(iii) Compute the Jacobian 7'"^'^'' = dLqjdq, where Lq denotes the set of equa- 
tions in operator form. 

(iv) Construct a simplified matrix A (preconditi oner), which is eas y to invert, 
but still share the spectral properties of (IHackbuschl . 1 1 9941) . 



1.0 



o.t 



Explicit 





10000.0 


\ 
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10-* 



to- 



10= 



100 



1000 



Fig. 1.6. The profiles of the CFL- number (solid line) versus the number of iteration both 
for explicit and implicit solution procedures (dashed line). The profiles correspond to the 
free-fall of spherical plasma onto a non-magnetized Schwarzschild black hole, in which 
the final solution is time-independent. 



(v) Solve the system of equation: 



An - d. 



(1.17) 



where // is a vector of small correction, so that q^'^^ = q' + fi. 
In general A + 7™"', which implies that Equation 11.171 should be solved 
iteratively to assure that the maximum norm of the defect, ||J||cx), is suffi- 
ciently small. 
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We note that for sufficiently small 6t, the matrix II6t can be made similar to 
jreai^ hcncc they share the same spectral space. As a consequence, a variety of so- 
lution procedures can be constructed that range from purely explicit up to strongly 
implicit, depending on how similar the preconditioner A is to the real Jacobian. 
This naturally suggests the hierarchical solution scenario as a highly powerful nu- 
merical algorithm for enhancing t he robustness o f explicit schemes and optimizing 
their efficiency (Fig. 11.31 see also iHujeirali . 120051) 



1.5 Summary-I 

In this part of the review we have presented a method for converting conditionally- 
stable explicit methods into numerically stable implicit solution procedures. The 
conversion method allows a considerable enlargement of the range of application 
of explicit methods. The hierarchical solution scenario is best suited for gradual 
enhancement of their robustness and optimizing their efficiency. 



1.6 Why Boltzmann Solvers? 



17 



Part II 

(Magneto-)Hydrodynamic Boltzmann Solvers 

In this part, I will discuss the implementation of a time-explicit gas-kinetic grid- 
based integrator for non-relativistic hydrodynamics introduced by Prendergast & 
Xu (1993), Xu (1999) and Tang & Xu (2000), and its extension to non-ideal 
magneto-hydrodynamics (Heitsch et al. 2004, 2007). Some properties of Boltz- 
mann solvers are discussed in ^1.6[ the equations and the implementation are de- 
scribed in 31.71 followed by a selection of test cases and applications ( 31.81 ) and a 
summary ( 31.91 ). 



1.6 Why Boltzmann Solvers? 

It is the physical model for the fluid equations which distinguishes gas-kinetic 
schemes from the widely popular Godunov methods. The latter are formulated 
on the basis of the Vlasov-equation, i.e. assuming that any dynamical time scale 
is lai^ger than the collision time between particles, setting the collision term in 
the Boltzmann equation to zero. The distribution function is then given by a 
Maxwellian at all times. In contrast, gas-kinetic schemes keep the collision term 
in the Boltzmann equation, but because of the impractibility to compute all the 
collisions between particles, they need to come up with a model for the collision 
term. 

One such model has been introduced by Bhatnagar, Gross & Krook (1954), for- 
mulating the collision term as the difference between the equilibrium distribution 
function g (the Maxwellian) and the initial distribution function /, resulting in a 
Boltzmann equation of the form 

dtf + udj + uduf=^^, (1.18) 

T 

where t is the collision time. Integrating eq. ll.l8l over a time t gives (at position x) 

1 r' 

f{x,t,u) = - g{x-u{t-t'),t',u)e'^'~'^'''dt' +e-"^ fo{x-ut,0,u), (1.19) 
T Jo 

where t is the collision time, and /o the initial distribution function. For a complete 
description, see Xu (2001). Thus, the distribution function / at time t gets two 
contributions: one from the decaying initial conditions f{t = 0), and one from the 
growing equilibrium distribution g. 

The 0th, 1st and 2nd order velocity moments of the distribution function (here 
for a monatomic gas) 

8^p[-j exp(i(u-U)2) (1.20) 
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result in the (macroscopic) conserved quantities density p, momentum density pU 
and total energy density pE. The quantity A = m/(2kT). The corresponding mo- 
ments of the Boltzmann equation 11.181 give the conservation equations. The BGK 
collision term in eq. 1 1.1 81 gives then rise to a viscous flux, depending on the ratio 
of the CFL time step and a specified collision time. Thus, the Reynolds number of 
the flow can be controlled. The Prandtl number is 1 by construction. The scheme 
is upwind and it satisfies the entropy condition (Prendergast & Xu 1993, Xu 2001). 
The fully controlled dissipative term come at (close) to no extra computational 
cost. Fragmentation of hydrodynamically unstable systems due to numerical noise 
thus can be suppressed. Specifically, gas-kinetic schemes can easily provide a vis- 
cosity independent of grid geometry, thus allowing e.g. the modeling of disks on a 
cartesian grid (see Slyz et al. 2002). 

In the following I will discuss a specific implementation of a gas-kinetic solver, 
namely Proteus (see Heitsch et al. 2007). 



1.7 Equations and Implementation: Proteus 

Proteus solves the equations of non-ideal magnetohydrodynamics, with an Ohmic 
resistivity Aq,, and a shear viscosity v. 



dtP + V- ip\) = 



dtpy + V 



66 



62 



dtpE + V 



62 

pE\ + {p + — )v ■ 
8;r 



(v 6)6 

4jt 



= V n 



(1.21) 
(1.22) 



= V • (V • n) + AnJ^ (1.23) 



dfB + V ■ (v6 - 6v) = AnV^B, 



(1.24) 



The mechanism how to split the fluxes at the cell walls is described in detail 
by Xu (1999) and will not be repeated here. Viscosity and resistivity are imple- 
mented as dissipative fluxes. They require spatially constant coefficients Aq, and v. 
Ambipolar drift is implemented in the two-fluid description, currently only for an 
isothermal equation of state, though. 

Higher-order time accuracy is achieved by a TVD Runge-Kutta time stepping 
(Shu & Osher 1988). For second-order spatial accuracy, a choice of reconstruction 
prescriptions is available. 

Proteus offers two gas-kinetic solvers, the one just described, and a one-step 
integrator at 2nd order in time and space for hydrodynamics. The latter has been 
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Fig. 1.7. Logarithm of the damping rate (eq. 111.251 ) of a linear Alfven wave against loga- 
rithm of the Ohmic resistivity for k = 1,2,4. The resolution is N = 64. Lines denote the 
analytical solution. 

discussed in detail by Slyz & Prendergast (1999) and Slyz et al. (2005), so that we 
refer the interested reader to those papers. 

1.8 Test Cases and Applications 

1.8.1 ID: Resistively damped Linear Alfven Wave 

This one-dimensional test checks the resistive flux implementation as well as the 
accuracy of teh overall scheme. A linear- Alfen wave under weak Ohmic dissipation 
is damped at a rate of 



where is the Ohmic resistivity, and k - Ittk/L is the wave number of the Alfven 
wave, with k a natural number. The strongly damped case, where the decay dom- 
inates the time evolution, is uninteresting for our application, since the Ohmic re- 
sistivity is mainly used to control numerical dissipation. Figure 11.71 shows the 
damping rate against Ohmic resistivity /in for /c = 1 , 2, 4 at a grid resolution of 
N = 64. The damping rate is derived by measuring the amplitude of the wave at 
each full wave period. 

From Figure [LTl it is clear that, as one diminishes the value of Aq, there comes 
a point when the numerical resistivity of the scheme becomes comparable to the 
physical one, causing the measured damping rate to flatten out and depart from the 
analytical solution. For k = 4 and Aq = 0.1, the wave decays too quickly to allow 




(1.25) 
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a reliable measurement, and the system enters the strongly damped branch of the 
dispersion relation. However, we emphasize that even at 16 cells per wave length 
the resistivity range available to Proteus spans nearly two orders of magnitude. 



1.8.2 ID: Linear Alfven Waves in Weakly Ionized Plasmas 

The dispersion relation for a linear Alfven wave in a weakly ionized plasma splits 
into two branches (Kulsrud & Pearce 1969): a strongly coupled branch, for which 
the ion Alfven frequency Uk = kBj yj47Tpi «; v,„ = yp„, the ion-neutral collision 
frequency, and a weakly coupled branch, for which oj^ » y,„ yjpi/pn- The strongly 
coupled case leads to a dispersion relation of 

..4\l/2 ,.2 



CO 



-'T^' (1-26) 

4Vi„ 2Vin 



with 6 = pi/rhon- Thus, the strongly coupled Alfven wave travels at the neutral 
Alfven speed cah = B/ ^J4np^ and is increasingly damped with decreasing collision 
frequency. The weakly coupled branch leads to 

Now, the wave travels at the ion Alfven speed, and damping is proportional to v,„. 
Since ca„ ^{p^Jpi, the speeds can be widely disparate. 

Figure 11.81 shows the real and imaginary part of the Alfven wave frequency in a 
weakly ionized plasma. For simplicity, we vary the collision coefficient jad and 
keep the densities constant. Wave speed (upper panel) and damping term (lower 
panel) are well reproduced. 



1.8.3 2D: Current Sheet 

This test is taken from Gardiner & Stone (2005). A square domain of extent < 
x,y < 7. and of constant density po = 1 and pressure /jq - 0. 1 is permeated by 
a magnetic field along the y direction such that Byi^.i < x < 1.5) = -1, and 
By - 1 elsewhere. The ratio of thermal over magnetic pressure is /3 = 0.2. This 
setup results in two magnetic null lines, which then are perturbed by velocities 
Vx = vo sin(;ry). Here, we use an adiabatic exponent of y = 5/3 and employ the 
conservative formulation of the scheme. Figure 11.91 summarizes the test results in 
the form of the magnetic energy density {B^} against time. Different line styles 
stand for resistivities, and the line thickness denotes the model resolution. We ran 
tests at = 128^, 256^ and 512^. All models ran up to f = 4 and farther except for 
the 5 12^ -model at Aq = 0. A finite resistivity helps stabilizing the code. 
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Fig. 1.8. Logarithm of the frequency (upper panel) and damping rate (lower panel) for 
the linear Alfven wave in a partially ionized plasma. For simplicity, we vary the collision 
coefficient jad instead of the density. 




12 3 4 

lime 



Fig. 1.9. Current sheet test. Magnetic energy density {B^} against time. A finite resistivity 
Aq helps stabilize the code. Line thickness stands for resolution, line style for resistivity. 



The evolution of the system follows that described by Gardiner & Stone (2005), 
including the merging of magnetic islands until there are two islands per magnetic 
null line left, located approximately at the velocity anti-nodes. For zero resistivity 
(solid lines), the magnetic energy decay depends strongly on the resolution. This 
effect is reduced by increasing Aq. For logA^ = -5 (dashed lines), the energy 
evolution follows pretty much the curves for Aq, = (solid lines), indicating in- 
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sufficient resolution. For log /In - -4, the two higher resolutions start to separate 
from the lower resolution run, while at log /In = -3, the two higher resolutions 
lead to indistinguishable curves (dash-3dot lines). 

1.8.4 2D: Advection of a Field Loop 

A cylindrical current distribution (i.e. a field loop) is advected diagonally across the 
simulation domain. Again, we follow the implementation presented by Gardiner 
& Stone (2005). Density and pressure are both initially uniform at po - 1 and 
Pq - \, and the fluid is described as an ideal gas with an adiabatic exponent of 
y - 5/3. The computational grid at a resolution of x Ny = 128 x 64 extends 
over -1.0 < X < 1.0 and -0.5 < y < 0.5. The field loop is initialized via the 
z-component of the vector potential A, = ao(R - r), where ao = 10^^, R = 0.3 and 
r = {x^ + j^)^^^. The loop is advected at an angle of 30 degrees with respect to 
the X-axis. Thus, two round trips in x correspond to one crossing in y. Figure 11.101 
shows the initial magnetic energy density with the magnetic field vectors over- 
plotted (top), and the distribution after two time-units measured in horizontal 
crossing times (bottom). The overall shape is preserved, although some artifacts 
are visible. These results concerning the shape are similar to those of Gai^diner & 
Stone (2005), specifically, Proteus preserves the circular field lines. This test uses 
An = 0. 

The time evolution of the magnetic energy density corresponding to Figure [TTTOl 
is shown in Figure 11.111 Diamonds stand for Proteus results, the energy decay 
observed by Gardiner & Stone (2005) with ATHENA is indicated by the solid line, 
following their analytical fit. The energies are normalized to 1 . Clearly, Proteus is 
somewhat more diffusive. 

In summary, these numerical test cases demonstrate that Proteus models dissi- 
pative MHD effects accurately. Furthermore, it can advect geometrically complex 
magnetic field patterns properly. 

1.9 Summary 

Gas-kinetic schemes provide a robust and physical mechanism to solve the equa- 
tions of magneto-hydrodynamics. Dissipative effects can be fully controlled. I 
discussed a specific implementation of a gas-kinetic solver - Proteus -, includ- 
ing resistivity and (two-fluid) ambipolar diffusion. Details of the implementation 
have been presented elsewhere (Tang & Xu 2000, Heitsch et al. 2004, 2007), and 
an application to shear flows in magnetized fluids will be discussed by Palotti et 
al. (2008). 
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Fig. 1.10. Field loop advection test: magnetic energy density B~ at f = (top) and att - 2 
corresponding to two horizontal crossing times (bottom), with over-plotted field vectors. 
The grid resolution is x Ny - 128 x 64. 




Fig. 1.11. Normalized magnetic energy density against time (in units of horizontal cross- 
ing time) with the same parameters as in Figure [TTTOl Diamonds stand for Proteus results, 
and the energy evolution as observed in ATHENA is shown by the soUd Une. 
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